;******************************************
; plot time serie of temperature anoamaly and HW frequency 
;********************************************
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"
;********************************************
begin

;set input
 dirPath0 = ".../NFood_replication/data/" ; change to local path of your directory
 dirPath1 = dirPath0+"/yield_data"

 ;target variables for analysis
 ;CROPPROD1C: "1-yr grain product C (gridcell)"
 ;CROPPROD1C_HARVEST:long_name = "crop harvest flux to 1-yr grain product pool (gridcell)"
 ;CROPPROD1C_HARV_CFT:long_name = "crop harvest flux to 1-yr grain product pool for each crop (pft)"
 var = "CROPPROD1C_HARV_CFT" ; (time, cft, lat, lon)
 nvar = dimsizes(var);

 ;cases and years
 casename =(/"IRCP45Clm50BgcCropGs_cplhist_1deg_CNTL","IRCP45Clm50BgcCropGs_cplhist_1deg_RCP85LU","IRCP85Clm50BgcCropGs_cplhist_1deg_CNTL","IRCP85Clm50BgcCropGs_cplhist_1deg_RCP45CO2", "IRCP85Clm50BgcCropGs_cplhist_1deg_SAI01","IRCP85Clm50BgcCropGs_cplhist_1deg_MSB01","IRCP85Clm50BgcCropGs_cplhist_1deg_CCT01"/)
 ;case = (/"RCP45", "RCP85", "RCP85_SAI", "RCP85_MSB", "RCP85_CCT"/);
 case = (/"RCP45","RCP45_85LU","RCP85","RCP85_45CO2","RCP85_SAI","RCP85_MSB","RCP85_CCT"/)
 ncase = dimsizes(case);
 syr = 2006;
 fyr = 2100;
 csyr = sprinti("%0.4i",syr);
 cfyr = sprinti("%0.4i",fyr);
 nyr = fyr - syr + 1;

;read dimension
 fdim = addfile(dirPath0+"/land_cover/RCP85.1deg.landarea.nc", "r");
 ;time = fdim->time
 ;cft  = fdim->cft
 lat = fdim->lat
 lon = fdim->lon
 nlat = dimsizes(lat);
 nlon = dimsizes(lon);
 dims = (/nlat,nlon/);


 ;Eight active crops in CLM5:
 ; 17/18=Temperate Corn/irrigated, 19/20=Spring Wheat/irrigated, 23/24=Temperate Soybean/irrigated, 41/42=Cotton/Cotton_irrigated,
 ; 61/62=Rice/Rice_irrigated, 67/68=Sugarcane/Sugarcane_irrigated, 75/76=Tropical Corn/irrigated, 77/78=Tropical Soybean/irrigated
 ;Winter wheat (21/22) has not been implemented yet. Only Spring/summer wheat
 ;cft_c3_crop = 1; cft_temperate_corn = 3; cft_spring_wheat = 5 ... cft sequence in crop landunit
 pfts = (/17,18, 19,20, 23,24, 41,42, 61,62, 67,68, 75,76, 77,78 /)      ; crop index among all pfts
 cfts = (/17,18, 19,20, 23,24, 41,42, 61,62, 67,68, 75,76, 77,78 /)-15   ; cft data index
 ncft = dimsizes(cfts);




;=============read data for all cases
 do kk=0,ncase-1

  ;readin specific crop area weight
   crop_wgt = new((/nyr,ncft-4,nlat,nlon/),"float")
   crop_wgt0  = fbindirread(dirPath1+"/"+case(kk)+"_crop_area_ha_nyr.ncft.lat.lon", 0, (/nyr,ncft,nlat,nlon/), "float" ) ; 4d binary data
   crop_wgt   = crop_wgt0(:,0:11,:,:)
   crop_wgt(:,0,:,:) = crop_wgt0(:,0,:,:) + crop_wgt0(:,12,:,:) ;rainfed corn (temperate+tropical)
   crop_wgt(:,1,:,:) = crop_wgt0(:,1,:,:) + crop_wgt0(:,13,:,:) ;irrigated corn (temperate+tropical)
   crop_wgt(:,4,:,:) = crop_wgt0(:,4,:,:) + crop_wgt0(:,14,:,:) ;rainfed soybean (temperate+tropical)
   crop_wgt(:,5,:,:) = crop_wgt0(:,5,:,:) + crop_wgt0(:,15,:,:) ;irrig soybean (temperate+tropical)


   ;read yield data from binary for each case
   ;area sum yield for each crop: MtC/year [nyr,ncft,nlat,nlon]
   filename0 = dirPath1+"/"+case(kk)+"_yield_sum_MtC_nyr.ncft.lat.lon"
   yield_sum = fbindirread(filename0, 0, (/nyr,ncft,nlat,nlon/), "float")

   ;area average yield for each crop: tonnes/ha/year [nyr,ncft,nlat,nlon]
   filename1 = dirPath1+"/"+case(kk)+"_yield_avg_tC-ha_nyr.ncft.lat.lon"
   yield_avg = fbindirread(filename1, 0, (/nyr,ncft,nlat,nlon/), "float")

   ;cultivation area for each crop: [nyr,ncft,nlat,nlon] 
   filename2 = dirPath1+"/"+case(kk)+"_crop_area_ha_nyr.ncft.lat.lon"
   cft_area  = fbindirread(filename2, 0, (/nyr,ncft,nlat,nlon/), "float")


 
 ;==============save lat/lon data to ncdf files
   ;---Open a new NetCDF file to write to
   fout_name =  dirPath1+"/"+case(kk)+"_crop_data.2006-2100.nc"
   system("rm -f " + fout_name)
   fout = addfile(fout_name,"c")
 
   ;---Prepare input data variables with coordinate arrays attached
   ntim =  nyr  ; 2006-2100
   ncft =  ncft ; irrig/rainfed
   nlat =  nlat
   nlon =  nlon
   nvars = 3    ; production/yield/area
 
;   var_names = "var" + sprinti("%02i",ispan(1,nvars,1))
   time      =  todouble(ispan(syr,fyr,1))
   time!0    = "time" ;named variables
   time@long_name = "simulation start:end years"
   ;time@_FillValue = 1.e+36; default_fillvalue("double")
   cft       = todouble(pfts) 
   cft!0     = "cft"  ;named variables, same as name
   cft@long_name = "cft indices of 16 active crop types"
   ;time      = time
   lat       = lat
   lon       = lon 
   var_names = (/ "CROP_PROD","CROP_YIELD","CROP_AREA" /) 
   var_units = (/ "MtC/year","tonne C/ha/year","ha" /)
   var_type  = "double"
 
   ;---Define the dimension names and their sizes on the file
   dims          := (/ntim,ncft,nlat,nlon/)
   dim_names     = (/time!0,cft!0,lat!0,lon!0/)
   dim_unlimited = (/True,False,False,False/)
   filedimdef(fout,dim_names,dims,dim_unlimited)
 
   ;---Define each variable, its name, type, its dimension names 
   filevardef(fout,time!0,typeof(time),time!0)
   filevardef(fout,cft!0,typeof(cft),cft!0)
   filevardef(fout,lat!0,typeof(lat),lat!0)
   filevardef(fout,lon!0,typeof(lon),lon!0)
   do nv=0,nvars-1
     filevardef(fout,var_names(nv),var_type,dim_names)
   end do
 
   ;---Define each variable's attributes.
   filevarattdef(fout,time!0,time) ;Copies attributes from an input variable to one or more variables of a file
   filevarattdef(fout,cft!0, cft)
   filevarattdef(fout,lat!0, lat)
   filevarattdef(fout,lon!0, lon)

   ;make sure all variables have the same _FillValue
   ;otherwise error: NetCDF: Not a valid data type or _FillValue type mismatch
 
;   fatt = True ;logic variable has _FillValue = 1, cause mismatch error!
   fatt = lon  ;copy _FillValue from lon and then modify long_name and units
   do nv=0,nvars-1
     fatt@long_name = var_names(nv)
     fatt@units     = "units_"+var_units(nv)
     ;fatt@_FillValue = 1.e+36    ; not effective, it's preset by data type of fatt 
     ;fatt@missing_value = 1.e+36 ; default_fillvalue("double")
     filevarattdef(fout,var_names(nv),fatt)
   end do

   ;---NOW write the input variables to the output file.
   fout->time  = (/time/)
   fout->cft   = (/cft/)
   fout->lat   = (/lat/)
   fout->lon   = (/lon/)
   fout->$var_names(0)$ = (/yield_sum/)  ;production [MtC]
   fout->$var_names(1)$ = todouble((/yield_avg/))  ;yield [tC/ha]
   fout->$var_names(2)$ = todouble((/cft_area/))   ;crop area [ha]

   ;x = create_dummy_var(var_names(0),time,cft,lat,lon,var_type)
   ;do nv=0,nvars-1
   ;  fout->$var_names(nv)$ = (/x/)
   ;end do
 
   ;---Close file. Not necessary, but a good idea.
   delete(fout)



 end do ;ncase



;==============save summary data to csv files
;*************************************
 ;Global Total Annual Yield of each crop: MtC/year
 filename3 = dirPath1+"/"+"YieldTot_cft_case_MtC_year" ; [ncft,ncase,nyr]
 totyd_sum = fbindirread(filename3, 0, (/ncft,ncase,nyr/), "float")
 totyd_sum@_FillValue = default_fillvalue("float")

 ;Global average yield flux of each crop: tonne/ha/year
 filename4 = dirPath1+"/"+"YieldAvg_cft_case_tonne_ha_year"; [ncft,ncase,nyr]
 totyd_avg = fbindirread(filename4, 0, (/ncft,ncase,nyr/), "float" )
 totyd_avg@_FillValue = default_fillvalue("float")

 filename6 = dirPath1+"/"+"YieldAvg_60lat_cft_case_tonne_ha_year"
 totyd_avg_60lat = fbindirread(filename6, 0, (/ncft,ncase,nyr/), "float" )
 totyd_avg_60lat@_FillValue = default_fillvalue("float")

 ;Global crop area per cft per case for each year: Million hectare
 filename5 = dirPath1+"/"+"TotArea_cft_case_year_Mha"
 tot_area_cft = fbindirread(filename5, 0, (/ncft,ncase,nyr/), "float" )
 tot_area_cft@_FillValue = default_fillvalue("float")



 pfts = (/17,18, 19,20, 23,24, 41,42, 61,62, 67,68, 75,76, 77,78 /) ; 8 active crops, non-irrigated and irrigated
 cftname = (/ "Temperate Corn",  "Irrigated Temperate Corn", "Spring Wheat", "Irrigated Spring Wheat",  "Temperate Soybean", "Irrigated Temperate Soybean", "Cotton", "Irrigated Cotton", \
              "Rice", "Irrigated Rice", "Sugarcane", "Irrigated Sugarcane", "Tropical Corn", "Irrigated Tropical Corn", "Tropical Soybean", "Irrigated Tropical Soybean" /)


 ;reorganize crops for ploting
 cftname12 = (/ "Corn", "Irrigated Corn", "Wheat", "Irrigated Wheat", "Soybean", "Irrigated Soybean",  "Cotton", "Irrigated Cotton", "Rice", "Irrigated Rice", "Sugarcane", "Irrigated Sugarcane" /)
 ncft2 = dimsizes(cftname12);

 ;irrigated vs. rainfed for each crop ;
 totyd_sum_12 = new((/12,ncase,nyr/),"float");  ;6x2 CPFs irrigated vs. rainfed
 totyd_avg_12 = new((/12,ncase,nyr/),"float");  ;6x2 CPFs 

 ;sum temperate + tropical for corn (0/1) and soybean (4/5)
 totyd_sum_12 = totyd_sum(0:11,:,:) ;excl. tropical corn and soybean, add below
 totyd_sum_12(0,:,:) = totyd_sum(0,:,:)+totyd_sum(12,:,:) ;Rainfed corn (tropical+temperate)
 totyd_sum_12(1,:,:) = totyd_sum(1,:,:)+totyd_sum(13,:,:) ;irrigated corn
 totyd_sum_12(4,:,:) = totyd_sum(4,:,:)+totyd_sum(14,:,:) ;Rainfed soybean
 totyd_sum_12(5,:,:) = totyd_sum(5,:,:)+totyd_sum(15,:,:) ;irrigated soybean

 ;average temperate + tropical
 totyd_avg_12 = totyd_avg(0:11,:,:) ;excl. tropical corn and soybean, add below
 totyd_avg_12(0,:,:) = (totyd_sum(0,:,:)+totyd_sum(12,:,:))/(tot_area_cft(0,:,:)+tot_area_cft(12,:,:)) ;Rainfed corn (tropical+temperate)
 totyd_avg_12(1,:,:) = (totyd_sum(1,:,:)+totyd_sum(13,:,:))/(tot_area_cft(1,:,:)+tot_area_cft(13,:,:)) ;irrigated corn
 totyd_avg_12(4,:,:) = (totyd_sum(4,:,:)+totyd_sum(14,:,:))/(tot_area_cft(4,:,:)+tot_area_cft(14,:,:)) ;Rainfed soybean
 totyd_avg_12(5,:,:) = (totyd_sum(5,:,:)+totyd_sum(15,:,:))/(tot_area_cft(5,:,:)+tot_area_cft(15,:,:)) ;irrigated soybean

 ;crop area
 tot_area_cft_12 = tot_area_cft(0:11,:,:) ;excl. tropical corn and soybean, add below
 tot_area_cft_12(0,:,:) = tot_area_cft(0,:,:)+tot_area_cft(12,:,:) ;Rainfed corn (tropical+temperate)
 tot_area_cft_12(1,:,:) = tot_area_cft(1,:,:)+tot_area_cft(13,:,:) ;irrigated corn
 tot_area_cft_12(4,:,:) = tot_area_cft(4,:,:)+tot_area_cft(14,:,:) ;Rainfed soybean
 tot_area_cft_12(5,:,:) = tot_area_cft(5,:,:)+tot_area_cft(15,:,:) ;irrigated soybean

;======merge rainfed and irrigated per crop
 cftname6 = (/ "Corn", "Wheat", "Soybean", "Cotton", "Rice", "Sugarcane" /)

 totyd_sum_6 = totyd_sum_12(ispan(0,10,2),:,:) + totyd_sum_12(ispan(1,11,2),:,:)
 tot_area_cft_6 = tot_area_cft_12(ispan(0,10,2),:,:) + tot_area_cft_12(ispan(1,11,2),:,:)
 totyd_avg_6 = totyd_sum_6/tot_area_cft_6


;Save data to csv

  csv_file1 = dirPath1+"/cropyield_tot_cft12.csv"
  csv_file2 = dirPath1+"/cropyield_avg_cft12.csv"
  system("rm -rf " + csv_file1)
  system("rm -rf " + csv_file2)
  csv_file11 = dirPath1+"/cropyield_tot_cft6.csv"
  csv_file22 = dirPath1+"/cropyield_avg_cft6.csv"
  system("rm -rf " + csv_file11)
  system("rm -rf " + csv_file22)


; Make the 1D dims the same length as the total length of 3D [ncft,ncase,nyr],
; so we can write them to CSV file along with data. 

; for all years 2006:2100
  dims0 = dimsizes(totyd_sum_12) ; [ncft-4,ncase,nyr]
  cftname1d  = ndtooned(conform_dims(dims0,cftname12,0))
  case1d = ndtooned(conform_dims(dims0,case, 1))
  time1d = ndtooned(conform_dims(dims0,ispan(2006,2100,1), 2))

  dims2 = dimsizes(totyd_sum_6) ; [6 cft,ncase,nyr]
  cftname1d2  = ndtooned(conform_dims(dims2,cftname6,0))
  case1d2 = ndtooned(conform_dims(dims2,case, 1))
  time1d2 = ndtooned(conform_dims(dims2,ispan(2006,2100,1), 2))
  

;---Write header to CSV file.
  
  field_names2 = (/"Crop name", "Case", "Year", "Yield [tonne/ha/year]" /) ; (ncft,ncase,nyr,data)
  field_names1 = (/"Crop name", "Case", "Year", "Yield [MtC/year]", "Area [Mha]"/) ; 
  header1 = [/str_join(field_names1,",")/]
  header2 = [/str_join(field_names2,",")/]
  write_table(csv_file1, "w", header1, "%s")
  write_table(csv_file2, "w", header2, "%s")
  write_table(csv_file11, "w", header1, "%s")
  write_table(csv_file22, "w", header2, "%s")


;Write crop data to csv file
 ;---convert 3D array to 1D for writing to CSV
 totyd_sum_1d = ndtooned(totyd_sum_12)
 totyd_avg_1d = ndtooned(totyd_avg_12)
 tot_area_cft_1d = ndtooned(tot_area_cft_12)

 totyd_sum_1d_6 = ndtooned(totyd_sum_6)
 totyd_avg_1d_6 = ndtooned(totyd_avg_6)
 tot_area_cft_1d_6 = ndtooned(tot_area_cft_6)
 
 alist1  = [/cftname1d,case1d,time1d, totyd_sum_1d, tot_area_cft_1d/]
 alist2  = [/cftname1d,case1d,time1d, totyd_avg_1d/]
 alist11  = [/cftname1d2,case1d2,time1d2, totyd_sum_1d_6, tot_area_cft_1d_6/]
 alist22  = [/cftname1d2,case1d2,time1d2, totyd_avg_1d_6/]

 ;format = "%g,%g,%g,%g,%g"
 format1 = "%s,%s,%4d,%g,%g"
 format2 = "%s,%s,%4d,%g"

 write_table(csv_file1, "a", alist1, format1)
 write_table(csv_file2, "a", alist2, format2)
 write_table(csv_file11, "a", alist11, format1)
 write_table(csv_file22, "a", alist22, format2)


end
